function [x,y] = VarStack(Z,nx,mx,myx,ny,hx,x_t,y_tN,N)

% Control variables
y(1:N+1,1:ny) = [Z(1:ny,1:N) y_tN ]';

% Endogenous state variables
x1 = [x_t(1:mx,1)        Z(ny+1:ny+mx,1:N);
      x_t(mx+1:mx+myx,1) Z(1:myx,1:N)      ]';

% 1 if we have off-diagonal elements in the exogenous shocks, else 0 for
% uncorrelated shocks
nx1         = mx+myx;
hxOffDiag   = any(any(hx(nx1+1:end,nx1+1:end)-diag(diag(hx(nx1+1:end,nx1+1:end)))));

% Exogenous state variables
nx2         = nx-nx1;
x2          = zeros(1+N,nx2);
if hxOffDiag == 0
    % AR(1) processes
    for i=1:nx2
        x2(1,i)     = x_t(nx1+i,1);
        x2(2:N+1,i) = hx(nx1+i,nx1+i).^(1:N).*x2(1,i);
    end
else
    % VAR(1) processes
    x2_transpose      = zeros(nx2,N+1);
    x2_transpose(:,1) = x_t(nx1+1:end,1);
    for jj=2:1+N
        x2_transpose(:,jj) = hx(nx1+1:end,nx1+1:end)*x2_transpose(:,jj-1);
    end
    x2 = x2_transpose';
end
x = [x1, x2];

end

